//
// Author 
// Xianguo Lu
// lu@physi.uni-heidelberg.de 
//

#ifndef MATHFIT_H
#define MATHFIT_H


typedef Double_t (*BFUNC)(const Double_t *, const Double_t *);

class MathFit
{
public:
  //LS fit, x, y filled, independet error
  static void InputData(const TH1D *hh, const Double_t threshold=0);
  static TString MnFit(BFUNC fn, const Int_t npar, Double_t *gpar, Double_t *gerr, const Double_t *step_size=0x0, const Double_t *lowlim=0x0, const Double_t *hilim=0x0, const Bool_t *pfix=0x0);

 private:
  static void FitFCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  
  static void CalcChi2();
  static TString Kernel(const Double_t *inpar=0x0, const Double_t *step_size=0x0, const Double_t *lowlim=0x0, const Double_t *hilim=0x0, const Bool_t *pfix=0x0);

  //______________________________________________
  //______________________________________________
    
  //------ Initialized in MnFit(), func, npar, par and err passed in from outside ------
  //------ all reset after Kernel ------
  static BFUNC fFunc;

  static Int_t fNPar; 
  static Double_t *fPar;//fPar[fNPar] = fChi2, fPar[fNPar+1]=fNDF
  static Double_t *fErr;
  
  //chi2, ndf only used internally, since if it stays, it may be confused by different mnfit. Reset after Kernel and only passed out by TString output
  static Double_t fChi2;
  static Int_t fNDF;

  static Int_t fNStep;
  
  //------ initialized via InputData() ------
  static Int_t fNData;
  
  static Double_t *fX;
  static Double_t *fY;
  
  static Double_t *fYerr;

  //------ Initialized via constructor -----
  const static Int_t fMaxNData;
  
  //------ initialized and destroied in Kernel
  static TMinuit * fMinuit;
  
  //------
};

//----------------------------------------------------------------

BFUNC MathFit::fFunc = 0x0;
Int_t MathFit::fNPar = 0;
Double_t * MathFit::fPar = new Double_t[20];
Double_t * MathFit::fErr = new Double_t[20];
Double_t MathFit::fChi2 = -999;
Int_t MathFit::fNDF = -999;
Int_t MathFit::fNStep = 0;

Int_t MathFit::fNData = 0;
const Int_t MathFit::fMaxNData = 1000;
Double_t * MathFit::fX = new double[1000];
Double_t * MathFit::fY = new double[1000];
Double_t * MathFit::fYerr = new double[1000];
TMinuit * MathFit::fMinuit = 0x0;

void MathFit::InputData(const TH1D *hh, const Double_t threshold)
{
  fNData = 0;

  fX[0] = -999;
  fY[0] = -999;
  fYerr[0]=-999;

  //------

  const TAxis * ax = hh->GetXaxis();
  for(Int_t id=ax->GetFirst(); id<=ax->GetLast(); id++){
    const Double_t yy=hh->GetBinContent(id);
    if(yy <= threshold)//"=" important to exclude 0-bins
      continue;

    fX[fNData] = hh->GetBinCenter(id);
    fY[fNData] = yy;

    fYerr[fNData]=hh->GetBinError(id);
   
    fNData++;
    if(fNData>=fMaxNData){
      printf("exitreport MathFit:: Ini fNData>=fMaxNData !! %d\n", fMaxNData);
      exit(1);
    }
  }
}

void MathFit::CalcChi2()
{
  fChi2 = 0;
  for(Int_t idata=0; idata<fNData; idata++){
    const Double_t evaly = fFunc(&(fX[idata]), fPar);
    const Double_t dy = fY[idata] - evaly;
    
    const Double_t den = pow(fYerr[idata],2);
    const Double_t part = dy*dy/den;
    
//    for(Int_t ip=0; ip<fNPar; ip++){
//      printf("par%d %e\n", ip, fPar[ip]);
//    }
//    printf("%d : %e %e %e -- %e %e -- %e %e %e -- %e\n", idata, fX[idata], fY[idata], fYerr[idata], evaly, dy, dfdx, den, part, fChi2);

    fChi2 += part;
  }
}

TString MathFit::Kernel(const Double_t *inpar, const Double_t *step_size, const Double_t *lowlim, const Double_t *hilim, const Bool_t *pfix)
{
  for(Int_t idata=0; idata<fNData; idata++){
    printf("%7d %13.6e %13.6e %10.2e\n", idata, fX[idata], fY[idata], fYerr[idata]);
  }
  
  printf("\n================= Starting MathFit ndata %d =================\n\n", fNData);
  
  fMinuit = new TMinuit(fNPar);
  fMinuit->SetFCN(FitFCN);
  
  const Double_t errdef = 1;//LS
  fMinuit->SetErrorDef(errdef);

  for(Int_t i=0;i<fNPar;i++){
    const Double_t start=inpar?inpar[i]:0.5;
    //here is only INITIAL step size
    const Double_t ss=step_size?step_size[i]:0.01;
    const Double_t low=lowlim?lowlim[i]:0;
    const Double_t hi=hilim?hilim[i]:0;

    fMinuit->DefineParameter(i, Form("p%d",i), start, ss, low,hi);

    if(pfix && pfix[i]){
      fMinuit->FixParameter(i);
    }
  }


  const Int_t maxiter = 5000;  
  //not really working
  //fMinuit->SetMaxIterations(maxiter);
  //fMinuit->Migrad();

  //tolerance = 0.1 by default; see http://www.dnp.fmph.uniba.sk/cernlib/asdoc/minuit/node18.html  (the MIGrad entry)
  const Double_t tol= 1;
  Double_t arg[2]={maxiter, errdef*tol};
  Int_t migradflag = -999;
  fMinuit->mnexcm("MIGRAD", arg, 2, migradflag);

  //improve error bars a lot! especially in ML fit!
  //only 1 argument: maxcalls
  Int_t hesseflag = -999;
  fMinuit->mnexcm("HESSE", arg ,1, hesseflag);

  //MINOS has asy. error
  
  //pass out par even not converged. Important because for pre-fit it may not converge, but the parameters are still useful
  for(Int_t ipar=0;ipar<fNPar;ipar++){
    fMinuit->GetParameter(ipar,fPar[ipar],fErr[ipar]);
  }
  fNDF = fNData-fMinuit->GetNumFreePars();
  fPar[fNPar]   = -999;
  fPar[fNPar+1] = fNDF;

  //4: abnormal termination (e.g., MIGRAD not converged)
  if(migradflag+hesseflag){
    printf("\nMathFit::Kernel fitting error !! %d %d\n\n", migradflag, hesseflag);
  }
    
  CalcChi2();
  printf("======> chisq/NDOF: %f/%d\n",fChi2,fNDF);
  
  fPar[fNPar] = fChi2;

  TString info(Form("%.3e %.0f err%d %d %d", fPar[fNPar], fPar[fNPar+1], migradflag+hesseflag, migradflag, hesseflag));

  delete fMinuit;
  fMinuit = 0x0;

  return info;
}

TString MathFit::MnFit(BFUNC fn, const Int_t npar, Double_t *gpar, Double_t *gerr, const Double_t *step_size, const Double_t *lowlim, const Double_t *hilim, const Bool_t *pfix)
{
  if(fNData==0){
    printf("MathFit::MnFit no data yet, do Ini!!\n");
    exit(1);
  }
  
  fFunc = fn;
  fNPar = npar;
  fPar=gpar;
  fErr=gerr;
  
  const Double_t *inpar=fPar;
  
  // Start the stopwatch. If reset is kTRUE reset the stopwatch before starting it (including the stopwatch counter). Use kFALSE to continue timing after a Stop() without resetting the stopwatch.
  TStopwatch stwa;
  stwa.Start(kTRUE);
  TString info =  Kernel(inpar , step_size, lowlim, hilim, pfix);
  stwa.Stop();
  stwa.Print("m");//"m": milli sec; "u": micro sec

  fFunc = 0x0;
  fNPar = 0;
  fPar = 0x0;
  fErr = 0x0;

  fNStep=0;

  fChi2=0;
  fNDF=0;

  return info;
}

//--------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------

void MathFit::FitFCN(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
  fNStep++;

  for(Int_t ii=0; ii<fNPar; ii++){
    fPar[ii] = par[ii];
  }

  CalcChi2();
  f= fChi2;

//  printf("%d %f %d ---",ptr->fNStep, f, npar);
//  for(Int_t idp=0; idp<npar; idp++)
//    printf("%15.6e ", par[idp]);
//  printf("\n");

}

#endif
